Density functional theory of dissipative systems 
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Much recent interest has focused on using single 
molecules as transistors for a new breed of computersyj. 
The complex nature of these devices, especially the leads, 
suggests that their properties can be sensitive to chemical 
details. Thus we wish to model the transport character- 
istics of such devices with first-principles electronic struc- 
ture methods, such as density functional theory (DFT). 
However, the traditional theorems of DFT do not apply 
to systems carrying current in finite electric fields. 

Valid applications of DFT are derived from exact prin- 
ciples of quantum mechanics. These involve proofs of a 
one-to-one correspondence between densities and poten- 
tials, i.e., that a given one-electron density can only be 
produced by at most a single one-body potential, under 
a given set of restrictions. Ground-state DFT was estab- 
lished by Hohenberg and Kohn(3|, by proving that, for 
interacting electrons, a given ground-state density can be 
produced by at most one ground-state one-body poten- 
tial. We say the one-body potential, and hence all other 
properties, is a functional of the ground-state density. 

Similarly, Runge and Gross (RG) showed^ that, for 
interacting electrons in a given initial state, a given evolu- 
tion of the one-body density can be produced by at most 
one time-dependent one-body potential, thus establishing 
the validity of time-dependent density functional theory 
(TDDFT). We usually start in a non-degenerate ground- 
state, so that the initial-state dependence becomes a 
ground-state density dependence, via Hohenberg-Kohn. 
A slight variation on this, proven in passing by RG, uses 
the one-body current density as the basic variable, called 
time-dependent current-density functional theory (TD- 
CDFT). This yields a more natural description of solids 
in uniform electric fields^], including the polarization- 
dependence first identified by GGGj5j]. 

On the other hand, transport problems within tradi- 



tional wavefunction theory are usually handled within the 
Landauer formalism^. The molecule and contacts are 
placed between two infinite reservoirs at chemical poten- 
tials that differ by the voltage drop across the molecule. 
A standard integral over the Green's function and cou- 
pling to the reservoirs, and Fermi occupation factors, 
then yields the current, and thus the conductance. This 
can be made exact by using non-equilibrium Keldysh 
Green's functions^, by imagining the system begins in 
the distant past with reservoirs and molecule decoupled, 
and the coupling turned on adiabatically to produce a 
steady-state current. Unfortunately, as the decoupled 
initial state is not the ground state of any Hamiltonian, 
the functional needed in TDCDFT differs in some un- 
known way from the usual functional Q. 

The present state-of-the-art of DFT transport through 
single molecules is embodied by calculations such as those 
of Rcfs. ^ and |^, and many others. A self-consistent 
Kohn-Sham (KS) calculation is performed for a molecule 
trapped between two leads, and the scattering states 
are inserted into the (two-terminal) Landauer formula 
to calculate the conductance. At finite fields, the KS 
Green's function is significantly distorted from its zero- 
field value. As is well-recognized j6j, use of the ground- 
state KS Green's function in place of the exact non- 
equilibrium Green's function is an unjustified approxima- 
tion. In particular, when coupling between the molecule 
and leads is weak, the poles of this Green's function, 
representing resonances of the molecule, are at the KS 
orbital energy differences, which are not the true exci- 
tations of the system. Acknowledging such limitations, 
there have been several recent attempts to go beyond this 
picture[il[ill. 

In the present paper, we prove a new density func- 
tional theorem that encompasses transport at finite elec- 
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trie fields, by including dissipation to phonons via a mas- 
ter equation. We derive the associated KS master equa- 
tion. The new exchange-correlation (XC) potential re- 
duces to that of TDCDFT in the limit of zero dissipation. 
Finally, we show how all the ingredients for a realistic 
calculation can be constructed. 

We consider only symmetric leads. To avoid the use 
of reservoirs, put the entire system on a long thin ring, 
and thread through the center a solenoidal magnetic field. 
This produces a spatially uniform electric field through- 
out the entire system, and is equivalent to a change of 
gauge [l^. The system is finite, and nowhere are there 
two different chemical potentials. In the limit of zero 
bias, the Kubo response formalism can be applied, in ei- 
ther traditional many-body theory or DFT. Using this 
geometry, Kohn and Kamenev'lS^ showed how, by being 
careful with the order of limits, one can, within time- 
dependent Hartree theory, consider the response of the 
system to an AC electric field of frequency and by 
taking w — > 0, recover the two- and four-terminal Lan- 
dauer results. This procedure has been generalized to 

tdcdftQ. 

However, an important difficulty arises when the field 
is finite. For a purely electronic system, the electrons will 
accelerate indefinitely, and the current grow infinitely. 
In reality, there is dissipation due to scattering with 
phonons to bring the system to equilibrium. This ba- 
sic phenomenon is described by the quantum Liouville 
equation for the density matrix of the entire system of 
electrons and phonons, Stot- 



dt 



= -i[Htot, Stotit)] 



(1) 



In the case of bulk transport, Kohn and Luttinger[l5| 
showed how, for scattering from dilute impurities in weak 
fields, Eq. |^ recovers the Boltzman equation, identify- 
ing the diagonal elements of the electronic density matrix 
with the distribution function. 

In the quantum mechanics of dissipative systems, there 
is a well-established procedure for incorporating the ef- 
fects of inelastic scattering with a reservoir into the long- 
time evolution of the system. In this case, the total 
Hamiltonian consists of the system Hamiltonian H (the 
electrons) and a reservoir Hamiltonian R (the phonons), 
coupled by K. The system Hamiltonian contains N inter- 
acting electrons. The coupling is linear in the phonon- 
coordinates, and involves only one-body forces on the 
system. Before an initial time {t = 0), both electrons 
and phonons are in thermal equilibrium at temperature 
T. The exact system density matrix S then satisfies 



dS{t) 
dt 



-i[H{t),S{t)]-iTy[K,Stot{t)] 



(2) 



where the trace is over reservoir coordinates. To derive 
a master equation for S alone we coarse-grain over a 



time scale At that is long compared to electronic transi- 
tions and phonon correlations, but short compared to the 
relaxation time, i.e., the time-scale on which the electrons 
are losing energy to the reservoir. This yields 



dS{t) 
dt 



-i[H,~S{t)]+C[S{t)] 



(3) 



where C is a superoperator, found by applying Fermi's 
Golden rule to the scattering process, and determined 
by the coupling K and the reservoir spectral density. It 
is usually written in a basis of eigenstates of the time- 
independent (many-body) Hamiltonian, H\A) = Ea\A). 
In this basis, C[S] reads 



C[S] = - ^ ^A^B {LabLbaS+ 



A.B 



SLabLba — '^LbaSLai 



(4) 
(5) 



where the operators Lba represent a transition from 
state A to state -B, and the transition probabilities Ta^b 
are given by 



^a^b = 



'C>{lliab)\iab? [n{ujAB) 

'D{uJBA)\lAB?n[uJBA) 



1) Ea>Eb 
Ea < Eb 



(6) 

in terms of the clcctron-phonon coupling elements jab- 
'D(u)) is the density of states of phonons with frequency 
ijj [ijjAB = Ea — Eb), and n{uj) — 1/ (e'^^' — l) is the 
thermal occupation factor. The transition probabilities 
satisfy detailed balance 

e^X^{-EA/kBT)TA^B = eM-EB/kBT)TB^A, (7) 

so that the steady-state solution of Eq. ^ yields the 
thermal equilibrium density matrix of H . No matter 
what the initial density matrix, the steady-state solu- 
tion is always the same. Thus the master equation cou- 
ples statistical mechanics to quantum mechanics and, 
by including only secular contributions, goes beyond the 
Schrodinger equation for the electrons to build in irre- 
versible evolution. 

As long as the coupling is weak, the rate at which the 
system relaxes back to equilibrium (the relaxation time) 
will be much longer than all other scales, and the master 
equation applies. While there is much discussion about 
the validity of the master equation, and what physics 
is contained within it, in what follows we take Eq. Q 
as given, with all its merits and flaws, and show how to 
map the system to an effective single-particle system. Al- 
though the system Hamiltonian is time- independent, the 
time-dependence in the master equation is generated by, 
e.g., starting in a non-equilibrium density matrix, which 
evolves into the thermal equilibrium density matrix. This 
involves excitation and de-excitation of the electrons. 

Our goal is describe a many-electron system, evolv- 
ing under a master equation, by a collection of non- 
interacting electrons. To do this, we allow the one-body 
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potential of H to be time-dependent. This breaks the 
connection between C and H, since the final H might 
have an external electric field turned on, but the initial H 
not. We then construct a Runge-Gross style proof that, 
for fixed electron-electron interaction and superoperator 
C, and for a given initial density matrix 5*0, no two one- 
body potentials can give rise to the same time-dependent 
density n{rt). We assume that the potentials are Taylor- 
series expandable around t — Q, and that some coefficient 
in the expansion is not uniform in space. Begin with the 
equation of motion for the current density: 

^ ^ - 0«f ) 

= -tTr (j(r)[iJ,5])+Tr (j(r)C[5]) (8) 

where j(r) = ^ - S{r — rjpi -I- PiS{r — r.;) is the current- 
density operator. Using the fact that cyclic permutation 
does not alter the trace, we can write the first term on 
the right as iTr {[H,j{r)]S), which is the usual contri- 
bution from evolution under a Hamiltonian. Evaluating 
everything at t = 0, and considering two systems with 
possibly different potentials but the same initial density 
matrix and coupling, we find 



dA(j(r)) 



dt 



-no(r)AVvoxt(rt = 0) 



(9) 



just as in RG. Thus two potentials that differ at i = 
give rise to two different currents. 

One next shows that, for two systems whose initial 
Hamiltonians are the same, and two operators that are 
identical initially, but whose time-evolution differs. 



dA{A) 



dt 



= Trf^ 
t=o V dt 



(10) 



within the master equation, because both the commuta- 
tor and C are identical in both systems at t = 0. Ap- 
plying this result to the equation of motion for the k-th 
derivative of the currents, 



d^+'Ajirt) - -no(r) Vd^Av^^tirt) 



(11) 



where 9q = {d'' /dt'')\^_Q. Thus, any difference in any 
derivative of the potentials (other than a constant), pro- 
duces two different currents. This establishes a one-to- 
one correspondence between densities and currents. 

We could stop here, as it is the TD current DFT that is 
neededj^J for transport calculations. But, for generality, 
we also wish to establish a density functional theory, the 
final step of the RG theorem uses continuity to show that 
the densities must differ for two different potentials. This 
follows in the master equation evolution, since 



d(n(r)) 



dt 



t=o 



= -Tr (Vj(r)5(0)) +Tr 



(n(r)C[5(0)]) 

(12) 



and the last term is the same in both systems. The usual 
arguments about the vanishing of the potentials sufh- 
ciently rapidly at large distances then suffice for finite 
systems (l7^, or single- valuedness for periodic systems^. 
Breakdown of continuity in the master equation occurs 
because, in Eq. ifT^ , the superoperator provides a correc- 
tion to the usual statement under Hamiltonian evolution, 
i.e., some momentum is transferred to the reservoir. But 
use of Eq. ((T^ restores continuity, and the correction 
can even be written in terms of a current [l^ . 

We have established that the potential is a functional 
of the time-dependent current density for a given inter- 
action, statistics, initial density matrix, and coupling. In 
principle, we can apply the same argument with the inter- 
action set to zero, to produce a set of time-dependent KS 
equations whose one-body potential, Vg[n, Ss(0),C](rt), is 
defined to yield the exact n{rt), when evolved under the 
master equation. By subtracting the external potential 
and Hartree contribution, we find an XC potential that 
has the same dependencies as the KS potential, but also 
depends on the initial density matrix of the interacting 
system. Fortunately, we can subsume all dependence on 
the initial density matrices into the functional itself, by 
beginning in an equilibrium distribution for both the in- 
teracting and non-interacting systems. In that case, the 
Mermin functional, which is just a functional of the ini- 
tial density, determines the initial density matrix. 

Constructed this way, the KS system has certain 
pathologies. The superoperator in the many-body mas- 
ter equation is guaranteed to vanish only on the many- 
body density matrix, not on the KS density matrix. 
To compensate for this, the corresponding KS poten- 
tial might need to evolve forever, even after the KS sys- 
tem has settled into a steady state 0- Nevertheless, we 
can construct a practical KS scheme by constructing a 
KS superoperator, Cg. To do so, we define Vs{T){r) as 
the KS potential in the Mermin functional at tempera- 
ture T, i.e., the potential that, when thermally occupied 
with non-interacting electrons, reproduces the exact one- 
electron density at thermal equilibrium. We then apply 
perturbation theory for a weak interaction between non- 
interacting electrons in this potential and the phonons in 
the reservoir, yielding the analogs to Eqs. © above: 



_ . '^i^^j)\l^J? iniuJ^j) + 1) €, > 6^ 



where ei are the eigenvalues of — V^/2 -I- Ws(T)(r). The 
matrix elements are now evaluated for the interaction 
between the KS system and the bath. To find the KS 
master equation itself, we reduce the many-body Eq. Q 
to a single particle form by tracing out all other degrees 
of freedom, and using a Hartree-style approximation for 
the two-particle correlation functions appearing in C[S]. 
In the basis of the single-particle KS orbitals, \n),\m),\p), 
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we find[l£ 

dSn 



dt 



p 

P 

In this KS master equation, the steady-state equihbrium 
has a static potential. The approach rate, determined 
by 7, will not match that of the true system, but such 
effects are absorbed in the XC potential. The important 
point is that, if such a KS system exists, it is unique for 
the given coupling, by the theorem proven above. 

Next we discuss the XC functional, which depends on 
the coupling to the reservoir. One could imagine per- 
forming accurate wavefunction calculations for a uniform 
gas on a ring, with the given coupling, to produce a local 
density approximation for the master equation. But we 
argue that the usual approximations of TDDFT (or TD- 
CDFT, as needed), such as the adiabatic local density 
approximation, are likely to suffice, because the KS mas- 
ter equation includes dissipation, and drives the system 
to the thermodynamic steady state. Thus the effect of 
dissipation on XC is likely to be small, and might even 
vanish in the limit of weak coupling to the reservoir. The 
important XC effects are to correct the electronic tran- 
sition frequencies into the true transitions, and this is 
captured exactly by such an approximation. 

The dissipative part of the master equation, C[S] in 
Eqs. ^ and ((SJ depends on phonon frequencies uj, den- 
sities of state and coupling matrix elements 7„m- 
All these quantities can be extracted from first-principles 
density- functional linear response calculations 20] . This 
allows for a fully consistent DFT implementation of the 
dissipative dynamics. 

Lastly, we discuss the recovery of the Landauer result 
for the case of transport through a single molecule, in the 
limit of weak bias. This has recentlyT^I been derived 
using TDCDFT. Here we show the derivation of linear 
response for the more familiar DFT, and the current ver- 
sion is essentially the same (but more cumbersome). 

For a given KS master equation, assume the density 
matrix has evolved into its steady state, so that 



[Ho,So] = ~iCs{So) 



(15) 



If Cs has been constructed to thcrmalizc the eigenstates 
of Hq, both sides of this equation vanish. Now imagine 
perturbing the system with a weak time-dependent po- 
tential AV , which becomes constant after a finite time. 
Allow the system to relax back to its new steady state, 
with density matrix Sq + AS. Equating equal powers in 
the perturbation: 



where rj is the first derivative of Cs . Expanding all quan- 
tities in eigenstates of Hq and solving, we find 

IaVab — IbVab 



ASab = 



(17) 



Ea- Eb+ iriAB 

where 5*0 = i-e., /a is the Fermi occupa- 

tion factor of state A. Calculating the density change by 
tracing the density operator with AS', and recognizing 
that the system eigenstates are Slater determinants of 
orbitals, one recovers the usual density-density response 
function |23|. A similar derivation holds for the current- 
current response within TDCDFT. 

To understand what this means, consider a KS mas- 
ter equation with weak dissipation (rjij << LUij), on an 
infinite ring. Turn on a small but finite electric field, 
and evolve the system into a steady state. This is the 
isothermal conductivity found from the steady-state 
solution for the electrons coupled to the phonons. The 
derivation above shows that this reduces to the adia- 
batic conductivity as given by the Kubo response for- 
mula, which has recently been shown to recover the Lan- 
dauer result (with possible XC corrections! [l^ . 
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[Ho, AS] + tv[So] ■ AS = -[AV, So] 



(16) 
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